#ifndef GRID_TURBO
#define GRID_TURBO
#include <vector>
#include <fstream>
#include <map>
#include "math.h"

//Basic type definitions
//Vector type
class Vector {
public:
	double x;
	double y;
	double z;
	Vector(): x(0), y(0), z(0)		
	{
	};
	Vector(double _x, double _y, double _z) :
		x(_x),
		y(_y),
		z(_z)
		{
		};
	double mod()
	{
		return sqrt(x*x+y*y+z*z);
	};
	inline const Vector& operator+=(const Vector& a)
	{
		x += a.x;
		y += a.y;
		z += a.z;
		return *this;
	};
	inline const Vector& operator-=(const Vector& a)
	{
		x -= a.x;
		y -= a.y;
		z -= a.z;
		return *this;
	};
	inline const Vector& operator*=(const double& a)
	{
		x *= a;
		y *= a;
		z *= a;
		return *this;
	};	
	inline const Vector& operator/=(const double& a)
	{
		x /= a;
		y /= a;
		z /= a;
		return *this;
	};
};

inline double operator*(const Vector& a, const Vector& b)
{
	return (a.x*b.x + a.y*b.y + a.z*b.z);
};

inline Vector operator+(const Vector& a, const Vector& b)
{
	return Vector(a.x+b.x, a.y+b.y, a.z+b.z);
};

inline Vector operator-(const Vector& a, const Vector& b)
{
	return Vector(a.x-b.x, a.y-b.y, a.z-b.z);
};

inline Vector operator*(const double& a, const Vector& b)
{
	return Vector(a*b.x, a*b.y, a*b.z);
};

inline Vector operator*(const Vector& b, const double& a)
{
	return Vector(a*b.x, a*b.y, a*b.z);
};

inline Vector operator/(const Vector& b, const double& a)
{
	return Vector(b.x/a, b.y/a, b.z/a);
};

inline Vector operator-(const Vector& b)
{
	return Vector(-b.x, -b.y, -b.z);
};

inline Vector operator&(const Vector& a, const Vector &b)
{
	return Vector(a.y*b.z-b.y*a.z, b.x*a.z-a.x*b.z, a.x*b.y-a.y*b.x);
};



//Grid types
struct Node {
	int GlobalIndex;
	Vector P;	//Position
	int BMarker; //TO DO Delete
	int a;//useful parametr
};

struct Face {
	int GlobalIndex;
	int FaceCell_1;		//Index of cell 1
	int FaceCell_2;		//Index of cell 2
	Vector FaceCenter;	//Center of face
	Vector FaceNormal;	//Surface normal * square
	double FaceSquare;	//Face square
	std::vector<int> FaceNodes;	//Indexes of all face nodes
	int isExternal;		//If its an external face (FaceCell_2 == 0)
	int BCMarker;		//Marker of boundary condition type	
};

class Cell {
public:
	int GlobalIndex;		
	std::vector<int> Faces; //Indexes of all cell faces
	double CellVolume;		// Volume of cell
	Vector CellCenter;		// Center of cell
	double CellHSize;		// Size of cell	
	std::vector<int> Nodes; //Cell  node indexes
};


class Grid
{
public:
	//For ordered type of grid dimension sizes
	int nx;
	int ny;
	int nz;
	//And index map function (i,j,k) -> GlobalIndex
	std::map<std::vector<int>, int> ordered_index;	
	int GetIndexFromCoords(int i, int j, int k);
	std::map<int, std::vector<int> > global_index_to_coords;
	std::vector<int> GetCoordsFromIndex(int ind);	
	void AddNewOrderedInfo(Cell &c, int i, int j, int k);
	//Basic grid structure
	std::vector<Cell> cells;	//Array of cells
	std::vector<Face> faces;	//Array of faces
	std::vector<Node> nodes;	//Array of nodes
	//Boundary conditions
	//Incidency list ( local indexes for now)
	std::map<int, int> inc;
	//Map from Global Indexes to local
	std::map<int, int> face_index;
	std::map<int, int> cell_index;
	std::map<int, int> node_index;
	int MaxCellIndex;	//Maximal global index in cells array
	int MaxFaceIndex;	//Maximal global index in faces array
	int MaxNodeIndex;	//Maximal global index in nodes array
	void refresh_indexes()
	{
		MaxCellIndex = 0;
		MaxFaceIndex = 0;
		MaxNodeIndex = 0;
		cell_index.clear();
		for (int i = 0; i<cells.size(); i++) {
			cell_index[cells[i].GlobalIndex] = i;
			if (cells[i].GlobalIndex > MaxCellIndex) MaxCellIndex = cells[i].GlobalIndex;
		};
		face_index.clear();
		for (int i = 0; i<faces.size(); i++) {
			face_index[faces[i].GlobalIndex] = i;
			if (faces[i].GlobalIndex > MaxFaceIndex) MaxFaceIndex = faces[i].GlobalIndex;
		};
		node_index.clear();
		for (int i = 0; i<nodes.size(); i++) {
			node_index[nodes[i].GlobalIndex] = i;
			if (nodes[i].GlobalIndex > MaxNodeIndex) MaxNodeIndex = nodes[i].GlobalIndex;
		};
		return;
	};
	//Additional data about sum of normals if needed
	std::vector<Vector> normal;
	//Calculate all cell properties given Nodes
	void CalcCell(Cell& new_cell);
	//Calculate all face properties given Nodes and FaceCell_x
	void CalcFace(Face& new_face);
	//Clear grid data, free all possible memory
	void Clear();
	//Select subgrid from grid given bool vector ( is cell in subgrid or not)
	Grid SubGrid(std::vector<bool> isCell);	
	//Usefull functions
	void Deform(double factor);	//Deform grid
	void ConvertToCylinder(bool inlet, double shift, double inlet_size_x, double inlet_size_y);	//Adapt grid to represent cylinder
	void GridLoad(const char*);
	void GridLoadCyl(const char*);
	void GridLoad3DZhukov(const char *mfile, const char *trifile, const char *tetrfile, const char *coordfile);
	void ToTecPlot(const char*);
	void ToTecPlotHex(const char *fname);
	void GridLoad2D(char *, char *, char *);	
	void GridWrite(char *,char *, char *);
	void GridRead(const char *, const char *, const char *, int);
	void WriteInnerCells(char *,char *);
};

//Add new Cell.Global_Index <-> i,j,k connection
void Grid::AddNewOrderedInfo(Cell &c, int i, int j, int k) {
	std::vector<int> coords(3);
	coords[0] = i;
	coords[1] = j;
	coords[2] = k;
	ordered_index[coords] = c.GlobalIndex;
	global_index_to_coords[c.GlobalIndex] = coords;
	return;
}

//Get global index from coordinats
int Grid::GetIndexFromCoords(int i, int j, int k) {
	std::vector<int> coords(3);
	coords[0] = i;
	coords[1] = j;
	coords[2] = k;
	return ordered_index[coords];
}

//Get coordinates from global index
std::vector<int> Grid::GetCoordsFromIndex(int ind) {	
	return global_index_to_coords[ind];;
}

//Deform 2D grid by a given factor
void Grid::Deform(double factor)
{
	for (int i = 0; i<cells.size(); i++) {
		cells[i].CellVolume *= (factor * factor);
		cells[i].CellHSize *= factor;
		cells[i].CellCenter = factor * cells[i].CellCenter;
	};
	for (int i = 0; i<faces.size(); i++) {
		faces[i].FaceCenter = factor * faces[i].FaceCenter;
		faces[i].FaceNormal = factor * faces[i].FaceNormal;
		faces[i].FaceSquare *= factor;
	};
	for (int i = 0; i<nodes.size(); i++) {
		nodes[i].P = factor * nodes[i].P;
	};
	return;
};

//Convert 2D rectangular grid to represent cylinder with axis OX
void Grid::ConvertToCylinder(bool inlet, double shift, double inlet_size_x, double inlet_size_y) 
{
	double r_min = 1e200, r_max = 0;
	double x_min = 1e200, x_max = 0;
	double Pi = 3.1415926535897;
	double alpha = 0.2;
	std::map<int, bool> to_del;

	for (int i = 0; i<faces.size(); i++) {
		Face& f = faces[i];
		to_del[faces[i].GlobalIndex] = false;
		Node& n1 = nodes[node_index[f.FaceNodes[0]]];
		Node& n2 = nodes[node_index[f.FaceNodes[1]]];
		//Vertical face
		if (n1.P.x == n2.P.x) {
			double r1 = n1.P.y;
			double r2 = n2.P.y;
			f.FaceNormal = ((r2 + r1) * sin(2*alpha) / 2.0) * f.FaceNormal;
			f.FaceSquare *= (r2 + r1) * sin(2*alpha) / 2.0;
		};
		//Horizontal face
		if (n1.P.y == n2.P.y) {			
			if (n1.P.y != 0) {
				f.FaceNormal = 2 * (n1.P.y) * sin(alpha) * f.FaceNormal;
				f.FaceSquare *= 2 * (n1.P.y) * sin(alpha);
			} else {
				//Delete this face
				to_del[faces[i].GlobalIndex] = true;
				f.FaceNormal = 0 * f.FaceNormal;
				f.FaceSquare *= 0;//(1e-10);
			};
		};
	};	
	
	double inlet_square = 0;
	for (int i = 0; i<cells.size(); i++) {
		r_min = 1e200;
		r_max = 0;
		x_min = 1e200;
		x_max = 0;
		Vector sum(0,0,0);
		Cell& c = cells[i];
		for (int j = 0; j<c.Faces.size(); j++) {			
			Face& f = faces[face_index[c.Faces[j]]];
			Node& n1 = nodes[node_index[f.FaceNodes[0]]];
			Node& n2 = nodes[node_index[f.FaceNodes[1]]];
			sum = sum + f.FaceNormal;
			//Vertical face
			if (n1.P.x == n2.P.x) {
				if (n1.P.x > x_max) x_max = n1.P.x;
				if (n1.P.x < x_min) x_min = n1.P.x;				
			};
			//Horizontal face
			if (n1.P.y == n2.P.y) {
				if (n1.P.y > r_max) r_max = n1.P.y;
				if (n1.P.y < r_min) r_min = n1.P.y;				
			};
		};
		//Add 2 new faces
		Face f1, f2;
		f1.FaceSquare = c.CellVolume;
		f2.FaceSquare = c.CellVolume;
		f1.GlobalIndex = ++MaxFaceIndex;
		f2.GlobalIndex = ++MaxFaceIndex;
		f1.FaceNormal = Vector(0, -sin(alpha), -cos(alpha));
		f1.FaceNormal = f1.FaceSquare * f1.FaceNormal;
		f2.FaceNormal = Vector(0, sin(-alpha), cos(alpha));
		f2.FaceNormal = f2.FaceSquare * f2.FaceNormal;
		f1.BCMarker = 5;
		f2.BCMarker = 6;
		f1.FaceCell_1 = c.GlobalIndex;
		f1.FaceCell_2 = -1;
		f2.FaceCell_1 = c.GlobalIndex;
		f2.FaceCell_2 = -1;
		f1.isExternal = true;
		f2.isExternal = true;
		faces.push_back(f1);
		faces.push_back(f2);
		c.Faces.push_back(f1.GlobalIndex);
		c.Faces.push_back(f2.GlobalIndex);
		//Add one more face for inlet condition		
		if (inlet && (c.CellCenter.y >= inlet_size_y) && (c.CellCenter.x>=shift) && (c.CellCenter.x<=(shift + inlet_size_x))) {
			Face new_f, new_f2;
			new_f.FaceCell_1 = c.GlobalIndex;
			new_f.FaceCell_2 = -1;
			new_f.BCMarker = 7;
			new_f.FaceSquare = 2 * f1.FaceSquare / (Pi / alpha);
			inlet_square += new_f.FaceSquare;
			new_f.FaceNormal = new_f.FaceSquare * Vector(0,0,-1);
			new_f.FaceCenter = c.CellCenter;
			new_f.GlobalIndex = ++MaxFaceIndex;
			new_f.isExternal = true;
			new_f2 = new_f;
			new_f2.GlobalIndex = ++MaxFaceIndex;
			new_f2.FaceNormal = -new_f2.FaceNormal;
			new_f2.BCMarker = 8;
			new_f2.isExternal = true;
			faces.push_back(new_f);
			faces.push_back(new_f2);
			c.Faces.push_back(new_f.GlobalIndex);
			c.Faces.push_back(new_f2.GlobalIndex);
		};
		for (int j = 0; j<c.Faces.size(); j++) {
			if (to_del[c.Faces[j]]) c.Faces.erase(c.Faces.begin() + j);
		};
		//Adjust cell volume
		sum = sum + f1.FaceNormal + f2.FaceNormal;
		c.CellVolume = (sin(2*alpha) / 2.0) * (x_max - x_min) * (r_max*r_max - r_min*r_min);
	};

	std::vector<Face> old_faces(faces);
	faces.clear();
	for (int i = 0; i<old_faces.size(); i++) if (!to_del[old_faces[i].GlobalIndex]) faces.push_back(old_faces[i]);
	refresh_indexes();	
	return;
};


struct FaceIndexes
{	
	int i,j,k;
	FaceIndexes(int _i, int _j, int _k)
	{
		if (_i<_j) {
			i = _i;
			_i = _j;
			_j = i;
		};
		if (_i<_k) {
			i = _i;
			_i = _k;
			_k = i;
		};
		if (_j<_k) {
			i = _j;
			_j = _k;
			_k = i;
		};
		i = _i;
		j = _j;
		k = _k;
	};
};

bool operator<(const FaceIndexes& a, const FaceIndexes& b)
{
	if (a.i!=b.i) return a.i<b.i;
	if (a.j!=b.j) return a.j<b.j;
	if (a.k!=b.k) return a.k<b.k;
	return false;
};

void Grid::Clear() {
	cells.clear();	//Array of cells
	faces.clear();	//Array of faces
	nodes.clear();	//Array of nodes	
	//Incidency list ( local indexes for now)
	inc.clear();
	//Map from Global Indexes to local
	face_index.clear();
	cell_index.clear();
	node_index.clear();
	return;
};

//Calculate volume of tetrahedron
double CalcVolume(Vector r1, Vector r2, Vector r3, Vector r4) {
	return fabs((r2-r1) * ((r3-r1) & (r4-r1))) / 6;
};


//Calculate all cell properties given Nodes
void Grid::CalcCell(Cell& new_cell)
{
	//Calculate cell center
	std::vector<Vector> n(new_cell.Nodes.size());
	new_cell.CellCenter = Vector(0,0,0);
	for (int i = 0; i<new_cell.Nodes.size(); i++) {
		n[i] = nodes[node_index[new_cell.Nodes[i]]].P;
		new_cell.CellCenter += n[i];
	};
	new_cell.CellCenter /= new_cell.Nodes.size();
	//TO DO Unify
	//Calculate cell volume
	new_cell.CellVolume = 0;	
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[1], n[2]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[3], n[2]);

	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[1], n[5]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[4], n[5]);

	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[3], n[7]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[0], n[4], n[7]);

	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[1], n[2]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[3], n[2]);

	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[1], n[5]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[4], n[5]);

	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[3], n[7]);
	new_cell.CellVolume += CalcVolume(new_cell.CellCenter, n[6], n[4], n[7]);
	return;
};
//Calculate all face properties given Nodes and FaceCell_x
void Grid::CalcFace(Face& new_face)
{
	//TO DO Unify
	//Calculate face normal
	Vector r1 = nodes[node_index[new_face.FaceNodes[1]]].P - nodes[node_index[new_face.FaceNodes[0]]].P;
	Vector r2 = nodes[node_index[new_face.FaceNodes[2]]].P - nodes[node_index[new_face.FaceNodes[0]]].P;
	Vector r3 = nodes[node_index[new_face.FaceNodes[3]]].P - nodes[node_index[new_face.FaceNodes[0]]].P;
	new_face.FaceNormal = (r1 & r2) + (r2 & r3);
	new_face.FaceNormal /= 2;
	//Calculate face center
	new_face.FaceCenter = Vector(0,0,0);
	for (int i = 0; i<new_face.FaceNodes.size(); i++) {
		new_face.FaceCenter += nodes[node_index[new_face.FaceNodes[i]]].P;
	};
	new_face.FaceCenter /= new_face.FaceNodes.size();
	//Direct normal right
	Vector A = cells[cell_index[new_face.FaceCell_1]].CellCenter - new_face.FaceCenter;
	if (new_face.FaceNormal*A>0) {
		new_face.FaceNormal *= -1;
	};
	new_face.FaceSquare = new_face.FaceNormal.mod();
	//Link face to corresponding cells
	cells[cell_index[new_face.FaceCell_1]].Faces.push_back(new_face.GlobalIndex);
	normal[cell_index[new_face.FaceCell_1]] += new_face.FaceNormal;
	if (!new_face.isExternal) {
		cells[cell_index[new_face.FaceCell_2]].Faces.push_back(new_face.GlobalIndex);
		normal[cell_index[new_face.FaceCell_2]] -= new_face.FaceNormal;
	};
	return;
};

void Grid::GridLoad(const char *fname)
{
	std::ifstream ifs(fname);	
	//Load nodes
	int N_node;
	ifs >> N_node;
	nodes.resize(N_node);
	for(int i = 0; i<N_node; i++)
	{
		nodes[i].GlobalIndex = i;
		ifs >> nodes[i].P.x;
		ifs >> nodes[i].P.y;
		ifs >> nodes[i].P.z;
	}
	//Load cells and define connectivity list
	std::map<FaceIndexes, std::pair<int,int> > f;
	int N_cell;
	ifs >> N_cell;
	cells.resize(N_cell);
	faces.resize(0);
	Cell cell;
	int a1,a2,a3,a4;	
	std::vector<FaceIndexes> fi;
	for(int i=0; i<N_cell; i++)
	{
		cells[i].GlobalIndex = i;
		ifs >> a1;
		ifs >> a2;
		ifs >> a3;
		ifs >> a4;
		int marker;
		ifs >> marker;
		a1--; a2--; a3--; a4--;		
		cells[i].Nodes.resize(4);
		cells[i].Nodes[0] = a1;
		cells[i].Nodes[1] = a2;
		cells[i].Nodes[2] = a3;
		cells[i].Nodes[3] = a4;
		cells[i].CellCenter = 0.25*(nodes[a1].P+nodes[a2].P+nodes[a3].P+nodes[a4].P);
		Vector a = (nodes[cells[i].Nodes[1]].P - nodes[cells[i].Nodes[0]].P);
		Vector b = (nodes[cells[i].Nodes[2]].P - nodes[cells[i].Nodes[0]].P);
		Vector c = (nodes[cells[i].Nodes[3]].P - nodes[cells[i].Nodes[0]].P);
		cells[i].CellVolume = fabs((a*(b&c))/6.0);
		fi.clear();
		fi.push_back(FaceIndexes(a1,a2,a3));
		fi.push_back(FaceIndexes(a1,a2,a4));
		fi.push_back(FaceIndexes(a4,a2,a3));
		fi.push_back(FaceIndexes(a1,a4,a3));		
		for (int j = 0; j<4; j++) {			
			std::map<FaceIndexes, std::pair<int,int> >::iterator it = f.find(fi[j]);
			if ( it != f.end())
			{
				(*it).second.second = i;
			} else {
				f[fi[j]].first = i;
				f[fi[j]].second = -1;
			};
		};		
	};
	//Load external faces and their markers
	std::map<FaceIndexes, int> extf;
	extf.clear();
	int N_boundary;
	ifs >> N_boundary;
	for (int i = 0; i<N_boundary; i++) {
		int ind1,ind2,ind3, bc_marker;
		ifs>>ind1>>ind2>>ind3>>bc_marker;
		extf[FaceIndexes(ind1-1,ind2-1,ind3-1)] = bc_marker;
	};
	//Create all faces
	int i = 0;
	inc.clear();
 	std::map<FaceIndexes, std::pair<int,int> >::iterator it;
	for (it = f.begin(); it!=f.end(); it++) {		
		FaceIndexes ind = (*it).first;
		std::pair<int,int> link = (*it).second;
		Face new_face;
		new_face.GlobalIndex = i++;
		new_face.FaceCell_1 = link.first;
		new_face.FaceCell_2 = link.second;
		new_face.FaceNodes.clear();
		new_face.FaceNodes.push_back(ind.i);
		new_face.FaceNodes.push_back(ind.j);
		new_face.FaceNodes.push_back(ind.k);
		//Calc surface normal and center position given all nodes
		new_face.FaceCenter = (1.0/3.0)*
			(nodes[new_face.FaceNodes[0]].P
			+ nodes[new_face.FaceNodes[1]].P
			+ nodes[new_face.FaceNodes[2]].P);
		Vector A = cells[new_face.FaceCell_1].CellCenter - new_face.FaceCenter;
		Vector S1 = (nodes[new_face.FaceNodes[1]].P-nodes[new_face.FaceNodes[0]].P);
		Vector S2 = (nodes[new_face.FaceNodes[2]].P-nodes[new_face.FaceNodes[0]].P);
		Vector S = (S1&S2);
		if (S*A>0) {
			S.x = -S.x;
			S.y = -S.y;
			S.z = -S.z;
		};
		new_face.FaceSquare = S.mod()/2;
		new_face.FaceNormal = 0.5*S;		
		//Boundary markers if there are some
		if (link.second != -1) {
			//Inner surface
			new_face.isExternal = false;
			inc[new_face.FaceCell_1] = new_face.FaceCell_2;
			inc[new_face.FaceCell_2] = new_face.FaceCell_1;
			cells[new_face.FaceCell_1].Faces.push_back(new_face.GlobalIndex);
			cells[new_face.FaceCell_2].Faces.push_back(new_face.GlobalIndex);
		} else {
			//Boundary face
			new_face.isExternal = true;
			cells[new_face.FaceCell_1].Faces.push_back(new_face.GlobalIndex);
			new_face.BCMarker = extf[ind];
		};
		faces.push_back(new_face);
	};
	f.clear();
	fi.clear();
	extf.clear();
	refresh_indexes();
	return;
};

void Grid::GridLoadCyl(const char* file) 
{
	
	//Struct for index operating
	struct Index {
		int n1;
		int n2;
		int n3;
		Index(int _n1, int _n2, int _n3) {
			n1 = _n1;
			n2 = _n2;
			n3 = _n3;
		};
		int GetIndex(int i, int j, int k) {
			i = (i+n1) % n1; //Cyclic index
			if (j == 0) {
				return (n1 * n2 + 1) * k;
			} else {
				return (n1 * n2 + 1) * k + n1 * (j - 1) + i + 1;
			};
		};
	};

	std::ifstream ifs(file);	
	int n1, n2, n3;
	int ni, nj, nk;
	ni = nj = nk = 0;
	ifs>>n1>>n2>>n3;
	//Set dimensions
	nx = n1;
	ny = n2;
	nz = n3;
	Index ind(n1,n2,n3);
	int N = (n1 * n2 + 1) * (n3 + 1);
	nodes.resize(N);	
	//Read nodes
	for (int i = 0; i<N; i++) {
		int tmp;
		ifs>>nodes[i].P.x>>nodes[i].P.y>>nodes[i].P.z>>tmp;
		nodes[i].GlobalIndex = i;		
	};
	//Generate cells
	int k_front = 80;
	int n_cells = 0;
	for (int j = 1; j<=n2; j++) {
		for (int i = 0; i<n1; i++) {
			for (int k = 0; k<=(n3-1); k++) {
				//if (k<k_front) continue;
				n_cells++;
				Cell new_cell;
				new_cell.GlobalIndex = ind.GetIndex(i,j,k);				
				//Fill ordered_index structure
				AddNewOrderedInfo(new_cell, i, j, k);
				new_cell.Nodes.push_back(ind.GetIndex(i,j-1,k));
				new_cell.Nodes.push_back(ind.GetIndex(i,j,k));
				new_cell.Nodes.push_back(ind.GetIndex(i,j,k+1));
				new_cell.Nodes.push_back(ind.GetIndex(i,j-1,k+1));
				new_cell.Nodes.push_back(ind.GetIndex(i+1,j-1,k));
				new_cell.Nodes.push_back(ind.GetIndex(i+1,j,k));
				new_cell.Nodes.push_back(ind.GetIndex(i+1,j,k+1));
				new_cell.Nodes.push_back(ind.GetIndex(i+1,j-1,k+1));				
				cells.push_back(new_cell);
			};
		};
	};	
	//Generate faces
	//i - layer
	for (int j = 1; j<=n2; j++) {
		for (int i = 0; i<n1; i++) {
			for (int k = 0; k<=(n3-1); k++) {
				ni++;
				Face new_face;
				new_face.GlobalIndex = ind.GetIndex(i,j,k);				
				new_face.isExternal = false;
				new_face.FaceCell_1 = ind.GetIndex(i-1,j,k);
				new_face.FaceCell_2 = ind.GetIndex(i,j,k);
				new_face.FaceNodes.push_back(ind.GetIndex(i,j-1,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i,j,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i,j,k+1));
				new_face.FaceNodes.push_back(ind.GetIndex(i,j-1,k+1));				
				faces.push_back(new_face);
			};
		};
	};
	//j - layer
	for (int j = 1; j<=n2; j++) {
		for (int i = 0; i<n1; i++) {
			for (int k = 0; k<=(n3-1); k++) {
				nj++;
				Face new_face;
				new_face.GlobalIndex = ind.GetIndex(i,j,k) + N;					
				//Outer surface
				if (j == n2) {
					new_face.isExternal = true;
					new_face.FaceCell_1 = ind.GetIndex(i,j,k);
					new_face.FaceCell_2 = -1;
					new_face.BCMarker = 2;
				};
				//Inner face
				if (j != n2) {
					new_face.isExternal = false;
					new_face.FaceCell_1 = ind.GetIndex(i,j+1,k);
					new_face.FaceCell_2 = ind.GetIndex(i,j,k);;
				};
				new_face.FaceNodes.push_back(ind.GetIndex(i,j,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i+1,j,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i+1,j,k+1));
				new_face.FaceNodes.push_back(ind.GetIndex(i,j,k+1));				
				faces.push_back(new_face);
			};
		};
	};
	//k - layer
	for (int j = 1; j<=n2; j++) {
		for (int i = 0; i<n1; i++) {
			for (int k = 0; k<=n3; k++) {
				nk++;
				Face new_face;
				new_face.GlobalIndex = ind.GetIndex(i,j,k) + 2*N;	
				//Back border
				if (k == 0) {
					new_face.isExternal = true;
					new_face.FaceCell_1 = ind.GetIndex(i,j,k);
					new_face.FaceCell_2 = -1;
					new_face.BCMarker = 3;					
				};
				//Front border
				if (k == n3) {
					new_face.isExternal = true;
					new_face.FaceCell_1 = ind.GetIndex(i,j,k-1);
					new_face.FaceCell_2 = -1;
					new_face.BCMarker = 4;
				};
				//Inner face
				if ((k != 0) && (k != n3)) {
					new_face.isExternal = false;
					new_face.FaceCell_1 = ind.GetIndex(i,j,k-1);
					new_face.FaceCell_2 = ind.GetIndex(i,j,k);;
				};
				new_face.FaceNodes.push_back(ind.GetIndex(i,j-1,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i+1,j-1,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i+1,j,k));
				new_face.FaceNodes.push_back(ind.GetIndex(i,j,k));				
				for (int l = 0; l < new_face.FaceNodes.size(); l++) {
					int nInd = new_face.FaceNodes[l];
					nodes[nInd].BMarker = new_face.BCMarker;
				};
				faces.push_back(new_face);				
			};
		};
	};	
	refresh_indexes();	

	//Calculate all face and cell parameters	
	normal.resize(cells.size());
	for (int i = 0; i<cells.size(); i++) normal[i] = Vector(0,0,0);
	for (int i = 0; i<cells.size(); i++) CalcCell(cells[i]);
	for (int i = 0; i<faces.size(); i++) CalcFace(faces[i]);
	/*
	std::ofstream ofs2("normal.dat");
	for (int i = 0; i<cells.size(); i++) ofs2<<normal[i].x<<" "<<normal[i].y<<" "<<normal[i].z<<"\n";
	ofs2.close();
	std::ofstream ofs3("bad.txt");
	for (int i = 0; i<cells.size(); i++) if (normal[i].mod() > 0.0001) {
		int res = cells[i].GlobalIndex;
		int _k = res / (n1*n2 + 1);
		res = res / (n1*n2 + 1);
		int _j = res / n1 + 1; 
		int _i = res % n1;
		ofs3<<_i<<" "<<_j<<" "<<_k<<"\n";									
	};
	ofs3.close();
	//Output grid
	ifs.close();
	std::ofstream ofs("grid.dat");	
	ofs << "VARIABLES= \"X\", \"Y\", \"Z";	
	ofs << "\"\n";
	ofs << "ZONE T=\"D\" ,N="<< nodes.size() <<", E=" << cells.size()
		<< ", DATAPACKING = POINT, ZONETYPE = FEBRICK\n";
		 //<< ", F=FEVector, ET=BRICK" <<"\n";
	for( int k=0; k<nodes.size(); k++) {
		ofs << nodes[k].P.x << "     " << nodes[k].P.y << "     " << nodes[k].P.z ;		
		ofs << "\n";
	};
	ofs << "\n";
	for( int k=0; k<cells.size(); k++) {
		ofs << cells[k].Nodes[0]+1 << "     " << cells[k].Nodes[1]+1 <<"     " 
			<< cells[k].Nodes[2]+1 <<"     " << cells[k].Nodes[3]+1 << " "
			<< cells[k].Nodes[4]+1 << "     " << cells[k].Nodes[5]+1 <<"     " 
			<< cells[k].Nodes[6]+1 <<"     " << cells[k].Nodes[7]+1 << "\n";
	 }

    ofs.close();
	*/
	return;
};

void Grid::GridLoad3DZhukov(const char *mfile, const char *trifile, const char *tetrfile, const char *coordfile)
{
	std::ifstream ifs(mfile);	
	//Load nodes
	int N_cell;
	int N_tri;
	int N_node;
	ifs >> N_node;
	ifs >> N_cell;
	ifs >> N_tri;
	nodes.resize(N_node);
	ifs.close();
	ifs.open(coordfile);
	for(int i = 0; i<N_node; i++)
	{
		nodes[i].GlobalIndex = i;
		ifs >> nodes[i].P.x;
		ifs >> nodes[i].P.y;
		ifs >> nodes[i].P.z;
	}
	ifs.close();
	ifs.open(tetrfile);
	//Load cells and define connectivity list
	std::map<FaceIndexes, std::pair<int,int> > f;
	cells.resize(N_cell);
	printf("Hoho");
	faces.resize(0);
	Cell cell;
	int a1,a2,a3,a4;	
	std::vector<FaceIndexes> fi;
	for(int i=0; i<N_cell; i++)
	{
		cells[i].GlobalIndex = i;
		ifs >> a1;
		ifs >> a2;
		ifs >> a3;
		ifs >> a4;
		cells[i].Nodes.resize(4);
		cells[i].Nodes[0] = a1;
		cells[i].Nodes[1] = a2;
		cells[i].Nodes[2] = a3;
		cells[i].Nodes[3] = a4;
		cells[i].CellCenter = 0.25*(nodes[a1].P+nodes[a2].P+nodes[a3].P+nodes[a4].P);
		Vector a = (nodes[cells[i].Nodes[1]].P - nodes[cells[i].Nodes[0]].P);
		Vector b = (nodes[cells[i].Nodes[2]].P - nodes[cells[i].Nodes[0]].P);
		Vector c = (nodes[cells[i].Nodes[3]].P - nodes[cells[i].Nodes[0]].P);
		cells[i].CellVolume = fabs((a*(b&c))/6.0);
		fi.clear();
		fi.push_back(FaceIndexes(a1,a2,a3));
		fi.push_back(FaceIndexes(a1,a2,a4));
		fi.push_back(FaceIndexes(a4,a2,a3));
		fi.push_back(FaceIndexes(a1,a4,a3));		
		for (int j = 0; j<4; j++) {			
			std::map<FaceIndexes, std::pair<int,int> >::iterator it = f.find(fi[j]);
			if ( it != f.end())
			{
				(*it).second.second = i;
			} else {
				f[fi[j]].first = i;
				f[fi[j]].second = -1;
			};
		};		
	};
	ifs.close();
	ifs.open(trifile);
	//Load external faces and their markers
	std::map<FaceIndexes, int> extf;
	extf.clear();
	int N_boundary;
	ifs >> N_boundary;
	for (int i = 0; i<N_boundary; i++) {
		int ind1,ind2,ind3, bc_marker;
		ifs>>ind1>>ind2>>ind3>>bc_marker;
		extf[FaceIndexes(ind1-1,ind2-1,ind3-1)] = bc_marker;
	};
	//Create all faces
	int i = 0;
	inc.clear();
 	std::map<FaceIndexes, std::pair<int,int> >::iterator it;
	for (it = f.begin(); it!=f.end(); it++) {		
		FaceIndexes ind = (*it).first;
		std::pair<int,int> link = (*it).second;
		Face new_face;
		new_face.GlobalIndex = i++;
		new_face.FaceCell_1 = link.first;
		new_face.FaceCell_2 = link.second;
		new_face.FaceNodes.clear();
		new_face.FaceNodes.push_back(ind.i);
		new_face.FaceNodes.push_back(ind.j);
		new_face.FaceNodes.push_back(ind.k);
		//Calc surface normal and center position given all nodes
		new_face.FaceCenter = (1.0/3.0)*
			(nodes[new_face.FaceNodes[0]].P
			+ nodes[new_face.FaceNodes[1]].P
			+ nodes[new_face.FaceNodes[2]].P);
		Vector A = cells[new_face.FaceCell_1].CellCenter - new_face.FaceCenter;
		Vector S1 = (nodes[new_face.FaceNodes[1]].P-nodes[new_face.FaceNodes[0]].P);
		Vector S2 = (nodes[new_face.FaceNodes[2]].P-nodes[new_face.FaceNodes[0]].P);
		Vector S = (S1&S2);
		if (S*A>0) {
			S.x = -S.x;
			S.y = -S.y;
			S.z = -S.z;
		};
		new_face.FaceSquare = S.mod()/2;
		new_face.FaceNormal = 0.5*S;		
		//Boundary markers if there are some
		if (link.second != -1) {
			//Inner surface
			new_face.isExternal = false;
			inc[new_face.FaceCell_1] = new_face.FaceCell_2;
			inc[new_face.FaceCell_2] = new_face.FaceCell_1;
			cells[new_face.FaceCell_1].Faces.push_back(new_face.GlobalIndex);
			cells[new_face.FaceCell_2].Faces.push_back(new_face.GlobalIndex);
		} else {
			//Boundary face
			new_face.isExternal = true;
			cells[new_face.FaceCell_1].Faces.push_back(new_face.GlobalIndex);
			new_face.BCMarker = extf[ind];
		};
		faces.push_back(new_face);
	};
	f.clear();
	fi.clear();
	extf.clear();
	refresh_indexes();
	return;
};


void Grid::ToTecPlot(const char *fname)
{
	std::ofstream ofs(fname);	
	ofs << "VARIABLES= \"X\", \"Y\", \"Z";	
	ofs << "\"\n";
	ofs << "ZONE T=\"D\" ,N="<< nodes.size() <<", E=" << cells.size()
		 << ", F=FEVector, ET=TETRAHEDRON" <<"\n";
	for( int k=0; k<nodes.size(); k++) {
		ofs << nodes[k].P.x << "     " << nodes[k].P.y << "     " << nodes[k].P.z ;		
		ofs << "\n";
	};
	ofs << "\n";
	for( int k=0; k<cells.size(); k++) {
		ofs << cells[k].Nodes[0]+1 << "     " << cells[k].Nodes[1]+1 <<"     " 
			<< cells[k].Nodes[2]+1 <<"     " << cells[k].Nodes[3]+1 << "\n";
	 }

    ofs.close();
	return;
};


//Get subgrid from a grid
Grid Grid::SubGrid(std::vector<bool> isCell) {
	if (isCell.size() != cells.size()) {
		throw 1;
	};
	std::map<int,bool> facesleft;
	std::map<int,bool> nodesleft;
	facesleft.clear();
	nodesleft.clear();
	Grid subgrid;
	subgrid.faces.clear();
	subgrid.cells.clear();
	subgrid.nodes.clear();
	for (int i = 0; i<isCell.size(); i++) if (isCell[i]) {
		subgrid.cells.push_back(cells[i]);
		for (int j = 0; j<cells[i].Faces.size(); j++) facesleft[cells[i].Faces[j]] = true;
		for (int j = 0; j<cells[i].Nodes.size(); j++) nodesleft[cells[i].Nodes[j]] = true;
	};
	std::map<int, bool>::iterator it;
	for (it = facesleft.begin(); it!=facesleft.end(); it++) {
		if ((*it).second) subgrid.faces.push_back(faces[face_index[(*it).first]]);			
	};
	for (it = nodesleft.begin(); it!=nodesleft.end(); it++) {
		if ((*it).second) subgrid.nodes.push_back(nodes[node_index[(*it).first]]);
	};			
	subgrid.refresh_indexes();
	return subgrid;
};

//Sega part

Vector CircleCenter(Vector& A, Vector& B, Vector& C)
{
	double a = sqrt((B-C)*(B-C));
	double b = sqrt((A-C)*(A-C));
	double c = sqrt((B-C)*(B-C));
	double p = 0.5*(a + b + c);
	double S = sqrt(p*(p - a)*(p - b)*(p - c));
	double alfa_a = ((A-B)*(A-C));
	alfa_a *=a*a/(8*S*S);
	double alfa_b = ((B-A)*(B-C));
	alfa_b *= b*b/(8*S*S);
	double alfa_c = ((C-A)*(C-B));
	alfa_c *= c*c/(8*S*S);
	return (alfa_a*A) + (alfa_b*B) + (alfa_c*C);
}
// Calc a unit normal
Vector Find_Normal(const Vector& r)
{
	Vector n;
	n.z = 0;
	if(r.x==0) {n.y = 0; n.x = 1;}
	else {n.y = r.x/sqrt(r.x*r.x + r.y*r.y);
	n.x = -(r.y*n.y)/r.x;}
	return n;
}

//struct for storage date
struct Date
{
	int Index;
	int date[3];
};

//Find two common nodes from tiangles a & b and calculate distance
void Find_common(int& com1, int& com2, std::vector<Node> a, std::vector<Node> b, double& l_com)
{
	int c=0;
	int i_1,i_2;
	for(int i=0; i<=2; i++)
	{
		for(int j=0;j<=2;j++)
		{
			if((a[i].GlobalIndex==b[j].GlobalIndex)&&(c==0)){com1=a[i].GlobalIndex; i_1 = i; c++;break;}
			if((a[i].GlobalIndex==b[j].GlobalIndex)&&(c==1)){com2=a[i].GlobalIndex; i_2 = i; c++;break;}
		}if(c==2) break;
	}
	l_com = (a[i_1].P - a[i_2].P)*(a[i_1].P - a[i_2].P);
	return;
}

//transpose two quantities
void swap(int& x, int& y)
{
    int m = x;
    x = y;
    y = m;
}

//Finding nextdoor triangles function
void Find_Neigh(int Tr_In, std::vector<Date> Neigh, int *neigh)
{
	for(int i = 0; i <=2; i++) neigh[i] = Neigh[Tr_In].date[i];
	return;
}



class Triangle
{public:
	int Tr_Index;
	std::vector<Node> tr_node;
	Vector Tr_Centr;
	void Tr_Calc(int, std::vector<Date>,  std::vector<Cell>);
};

//Calculating triangle
void Triangle::Tr_Calc(int Tr_In, std::vector<Date> ele, std::vector<Cell> cells)
{
	Tr_Index = Tr_In;
	tr_node.resize(3);
	//Find vertices and vertices coordinates of triangle
	for(int i=0; i<=2; i++)
	{
	tr_node[i].GlobalIndex = ele[Tr_In].date[i];
	tr_node[i].P = cells[tr_node[i].GlobalIndex].CellCenter;
	tr_node[i].BMarker = (int)cells[tr_node[i].GlobalIndex].CellHSize;
	}
	//Calc the Center
	Tr_Centr = CircleCenter(tr_node[0].P, tr_node[1].P, tr_node[2].P);
	return;
}

void Grid::GridLoad2D(char *fname_ele,char *fname_node, char *fname_neigh)
{
	double read;
	//Read Voronoi Cells
	std::ifstream ifs_v_cells(fname_node);
	int N_cell;
	ifs_v_cells >> N_cell >> read >> read >> read;
	cells.resize(N_cell + 1);
	for(int i = 1; i<= N_cell; i++)
	{
		ifs_v_cells >> cells[i].GlobalIndex >> cells[i].CellCenter.x >> cells[i].CellCenter.y >> read;
		cells[i].CellCenter.z = 0;
		cells[i].CellVolume = 0;
		if(read==1) cells[i].CellHSize = 1;//fictions Cell
		if(read==0) cells[i].CellHSize = 0;
		if(read==2) cells[i].CellHSize = 2;//fictions Cell
	}
	ifs_v_cells.close();
	//usefull vectors for downloading date
	std::vector<Date> Neigh;
	std::vector<Date> Tr_Ele;
	//Loading of Neigh
	std::ifstream ifs_neigh(fname_neigh);
	int N_v_node;
	ifs_neigh >> N_v_node >> read;
	Neigh.resize(N_v_node+1);
	for(int i = 1; i<=N_v_node; i++)
	{
		ifs_neigh >> Neigh[i].Index >> Neigh[i].date[0] >> Neigh[i].date[1] >> Neigh[i].date[2];
	}
	ifs_neigh.close();
	//Loading of Triangle Ele
	std::ifstream ifs_tr_ele(fname_ele);
	ifs_tr_ele >> read >> read >> read;
	Tr_Ele.resize(N_v_node+1);
	for(int i = 1; i<=N_v_node; i++)
	{
		ifs_tr_ele >> Tr_Ele[i].Index >> Tr_Ele[i].date[0] >> Tr_Ele[i].date[1] >> Tr_Ele[i].date[2];
	}
	ifs_tr_ele.close();
			
	//Calc voronoi nodes
	nodes.resize(N_v_node+1);
	faces.resize(1);
	int k = 1;//Index of Edges
	for(int i = 1; i<=N_v_node; i++)
	{
		printf("%ld\n",i);
		nodes[i].a = 0;
		nodes[i].GlobalIndex = i;
		Triangle tr_our;
		tr_our.Tr_Calc(nodes[i].GlobalIndex, Tr_Ele, cells);
		nodes[i].P = tr_our.Tr_Centr;
		int neigh[3];
		//Find neighs
		Find_Neigh(nodes[i].GlobalIndex, Neigh, neigh);
		for(int j=0; j<=2; j++)
		{
			if(neigh[j]==-1)continue;//neigh[j] doesn't exist
			if(nodes[neigh[j]].a==1)continue;//We've considered this triangle
			Triangle tr_neigh;
			tr_neigh.Tr_Calc(neigh[j], Tr_Ele, cells);
			//Calculating new face
			Face New_face;
			New_face.FaceNodes.resize(2);
			Node face_cell1, face_cell2;
			double l_com;
			Find_common(face_cell1.GlobalIndex, face_cell2.GlobalIndex, tr_neigh.tr_node, tr_our.tr_node, l_com);
			New_face.GlobalIndex = k;
			New_face.FaceNodes[0] = tr_our.Tr_Index;
			New_face.FaceNodes[1] = tr_neigh.Tr_Index;
			New_face.FaceCenter = 0.5*(tr_our.Tr_Centr + tr_neigh.Tr_Centr);
			New_face.FaceSquare = sqrt((tr_our.Tr_Centr - tr_neigh.Tr_Centr)*(tr_our.Tr_Centr - tr_neigh.Tr_Centr));
			New_face.FaceNormal = Find_Normal(tr_our.Tr_Centr + tr_neigh.Tr_Centr);
			New_face.FaceNormal = New_face.FaceSquare*New_face.FaceNormal;
			double s_i = 0.25*New_face.FaceSquare*l_com;

			//calc Boundarity
			if((cells[face_cell1.GlobalIndex].CellHSize==0)&(cells[face_cell2.GlobalIndex].CellHSize==0))
			{
				New_face.FaceCell_1 = face_cell1.GlobalIndex;
				New_face.FaceCell_2 = face_cell2.GlobalIndex;
				New_face.isExternal = 0;
				New_face.BCMarker = 0;
				cells[face_cell1.GlobalIndex].CellVolume += s_i;
				cells[face_cell1.GlobalIndex].CellVolume += s_i;
				nodes[i].BMarker = 0;
			}else
			{
				if(cells[face_cell2.GlobalIndex].CellHSize!=0) swap(face_cell2.GlobalIndex,face_cell1.GlobalIndex);
				New_face.FaceCell_1 = face_cell2.GlobalIndex;
				New_face.FaceCell_2 = 0;
				New_face.isExternal = 1;
				New_face.BCMarker = (int)cells[face_cell1.GlobalIndex].CellHSize;
				cells[face_cell2.GlobalIndex].CellVolume += s_i;
				nodes[i].BMarker = New_face.BCMarker;
			}
			faces.push_back(New_face);
			k++;
		}
		nodes[i].a = 1;
	}
	refresh_indexes();
	return;
};
		//TO DO [i]->[i-1]
void Grid::GridWrite(char *fname_node, char *fname_cell, char *fname_edge)
{
	std::ofstream ofs_node(fname_node);	
	ofs_node << (nodes.size() - 1) <<	"\n";
	for(int i = 1; i<= nodes.size() - 1; i++)
	{
		ofs_node << nodes[i].GlobalIndex << "\t" << nodes[i].P.x << "\t" << nodes[i].P.y << "\t" << nodes[i].P.z << "\t" << nodes[i].BMarker<<"\n";
	}
	ofs_node.close();

	std::ofstream ofs_cell(fname_cell);	
	ofs_cell << (cells.size() - 1) <<	"\n";
	for(int i = 1; i<= cells.size() - 1; i++)
	{
		ofs_cell << cells[i].GlobalIndex << "\t" << cells[i].CellCenter.x << "\t" << cells[i].CellCenter.y << "\t";
		ofs_cell << cells[i].CellCenter.z << "\t" << cells[i].CellVolume <<"\n";
	}
	ofs_cell.close();

	std::ofstream ofs_edge(fname_edge);	
	ofs_edge << (faces.size() - 1) <<	"\n";
	for(int i = 1; i<= faces.size() - 1; i++)
	{
		ofs_edge << faces[i].GlobalIndex << "\t" << faces[i].FaceCell_1 << "\t" << faces[i].FaceCell_2<< "\t";
		ofs_edge << faces[i].FaceCenter.x << "\t" << faces[i].FaceCenter.y << "\t" << faces[i].FaceCenter.z << "\t";
		ofs_edge << faces[i].FaceNodes[0]<< "\t" << faces[i].FaceNodes[1] << "\t" << faces[i].FaceNormal.x << "\t";
		ofs_edge << faces[i].FaceNormal.y << "\t" << faces[i].FaceNormal.z << "\t" << faces[i].FaceSquare << "\t";
		ofs_edge << faces[i].isExternal << faces[i].BCMarker << "\t" <<"\n";
	}
	ofs_edge.close();
	return;
};

void Grid::GridRead(const char *fname_node, const char *fname_cell, const char *fname_edge, int Nodes_size)
{
	int N_node;
	std::ifstream ofs_node(fname_node);	
	ofs_node >> N_node;
	nodes.resize(N_node + 1);
	for(int i = 1; i<=N_node; i++)
	{
		ofs_node >> nodes[i].GlobalIndex >> nodes[i].P.x >> nodes[i].P.y >> nodes[i].P.z >> nodes[i].BMarker;
	}
	ofs_node.close();

	int N_cell;
	std::ifstream ofs_cell(fname_cell);
	ofs_cell >> N_cell;
	cells.resize(N_cell + 1);
	for(int i = 1; i<=N_cell; i++)				//rewrite this part for any number of Nodes
	{
		ofs_cell >> cells[i].GlobalIndex >> cells[i].CellCenter.x >> cells[i].CellCenter.y;
		ofs_cell >> cells[i].CellCenter.z >> cells[i].CellVolume >> cells[i].CellHSize;
		cells[i].Nodes.resize(Nodes_size);
		for(int j=0;j<Nodes_size;j++) ofs_cell >> cells[i].Nodes[j];
	}
	ofs_cell.close();

	int N_face;
	std::ifstream ofs_edge(fname_edge);	
	ofs_edge >> N_face;
	faces.resize(N_face + 1);
	for(int i = 1; i<=N_face; i++)
	{
		faces[i].FaceNodes.resize(2);
		ofs_edge >> faces[i].GlobalIndex >> faces[i].FaceCell_1 >> faces[i].FaceCell_2;
		ofs_edge >> faces[i].FaceCenter.x >> faces[i].FaceCenter.y >> faces[i].FaceCenter.z;
		ofs_edge >> faces[i].FaceNodes[0] >> faces[i].FaceNodes[1] >> faces[i].FaceNormal.x;
		ofs_edge >> faces[i].FaceNormal.y >> faces[i].FaceNormal.z >> faces[i].FaceSquare;
		ofs_edge >> faces[i].isExternal >> faces[i].BCMarker;
	}
	ofs_edge.close();
	//TO DO 
	for (int i = 0; i<faces.size()-1; i++) {
		faces[i] = faces[i+1];		
	};
	for (int i = 0; i<cells.size()-1; i++) {
		cells[i] = cells[i+1];
	};
	for (int i = 0; i<nodes.size()-1; i++) {
		nodes[i] = nodes[i+1];
	};
	faces.pop_back();
	cells.pop_back();
	nodes.pop_back();
	std::vector<Cell> old_cells(cells);
	cells.clear();
	for (int i = 0; i<old_cells.size(); i++) if (old_cells[i].CellHSize == 0) {
		cells.push_back(old_cells[i]);
	};
	refresh_indexes();
	for (int i = 0; i<faces.size(); i++) {
		cells[cell_index[faces[i].FaceCell_1]].Faces.push_back(faces[i].GlobalIndex);
		if (!faces[i].isExternal) cells[cell_index[faces[i].FaceCell_2]].Faces.push_back(faces[i].GlobalIndex);
	};
	printf("Data successfully loaded\n");
	return;
};

void WriteInnerCells(char *ele, char *node, char *inner) {
	std::vector<int> goodNodes(0);
	std::map<int, int> markers;
	std::ifstream ifs(node);
	int Nnode;
	int read;
	int n;
	double x;
	double y;
	int marker;
	ifs>>Nnode >> read >> read >> read;
	for (int i = 0; i<Nnode; i++) {
		ifs>>n>>x>>y>>marker;
		markers[n] = marker;
		if (marker == 0) goodNodes.push_back(n);
	};
	ifs.close();
	ifs.open(ele);
	int Nele;
	int n1, n2, n3;
	std::vector<std::vector<int> > innerTriangles(0);	
	ifs>>Nele>>read>>read;
	for (int i = 0; i<Nele; i++) {
		ifs>>n>>n1>>n2>>n3;
		if ((markers[n1] == 0) && (markers[n2] == 0) && (markers[n3] == 0)) {
			std::vector<int> T(3);
			T[0] = n1;
			T[1] = n2;
			T[2] = n3;
			innerTriangles.push_back(T);
		};
	};
	ifs.close();
	std::ofstream ofs(inner);
	ofs<<innerTriangles.size()<<"\n";
	for (int i = 0; i<innerTriangles.size(); i++) {
		for (int j = 0; j<3; j++) ofs<<innerTriangles[i][j]<<" ";
		ofs<<"\n";
	};
	ofs.close();
	return;
};

//Output solution to tecplot
void Grid::ToTecPlotHex(const char *fname)
{	
	std::ofstream ofs(fname);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Z\", \"pnode\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << nodes.size() << ", E=" << cells.size() <<", F=FEBLOCK, ET=BRICK\n";	
	ofs<<"VARLOCATION = (NODAL, NODAL, NODAL, CELLCENTERED)\n";
	//X
	for (int i = 0; i<nodes.size(); i++) {
		ofs<<nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<nodes.size(); i++) {
		ofs<<nodes[i].P.y<<"\n";
	};
	//Z
	for (int i = 0; i<nodes.size(); i++) {
		ofs<<nodes[i].P.z<<"\n";
	};
	
	//pnode
	for (int i = 0; i<cells.size(); i++) {
		ofs<<floor(cells[i].CellHSize)<<"\n";
	};
	
	//Connectivity list
	for (int i = 0; i<cells.size(); i++) {	
		ofs<<cells[i].Nodes[0]<<" "
			<<cells[i].Nodes[1]<<" "
			<<cells[i].Nodes[2]<<" "
			<<cells[i].Nodes[3]<<" "
			<<cells[i].Nodes[4]<<" "
			<<cells[i].Nodes[5]<<" "
			<<cells[i].Nodes[6]<<" "
			<<cells[i].Nodes[7]<<"\n";
	};
	ofs.close();	
	return;
};


#endif
